Core-halo distribution in the Hamiltonian Mean-Field Model 
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We study a paradigmatic system with long-range interactions: the Hamiltonian Mean-Field Model 
(HMF) . It is shown that in the thermodynamic limit this model does not relax to the usual equilib- 
rium Maxwell-Boltzmann distribution. Instead, the final stationary state has a peculiar core-halo 
structure. In the thermodynamic limit, HMF is neither ergodic nor mixing. Nevertheless, we find 
that using dynamical properties of Hamiltonian systems, it is possible to quantitatively predict both 
the spin distribution and the velocity distribution functions in the final stationary state, without 
any adjustable parameters. We also show that HMF undergoes a non-equilibrium first-order phase 
transition between paramagnetic and ferromagnetic states. 
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Since the early work of Clausius, Boltzmann, and 
Gibbs it has been known that for particles interacting 
through short-range potentials, the final stationary state 
reached by a system corresponds to the thermodynamic 
equilibrium pj. Although no exact proof exists, in prac- 
tice it is found that non-integrable systems with fixed 
energy and number of particles (microcanonical ensem- 
ble) always relax to a unique stationary state which only 
depends on the global conserved quantities: energy, mo- 
mentum, and angular momentum. The equilibrium state 
does not depend on the specifics of the initial particle dis- 
tribution. The situation is very different for systems in 
which particles interact through long-ranged unscreened 
potentials. This is the case for gravitational systems and 
confined one component plasmas 0, [3] ■ For these sys- 
tems, in the thermodynamic limit, collision duration time 
diverges, and the thermodynamic equilibrium is never 
reached [1]. Instead, as time i — > oo, these systems be- 
come trapped in a stationary state characterized by a 
broken ergodicity • Unlike the thermodynamic equi- 
librium, the stationary state depends explicitly on the 
initial particle distribution. Over the last 50 years, there 
has been a great effort to predict the final stationary 
state without having to explicitly solve the iV-body dy- 
namics or the coUisionless Boltzmann (Vlasov) equation. 
Qualitatively, it has been observed that for many dif- 
ferent systems the non-equilibrium stationary state has 
a peculiar core-halo shape. Recently, an ansatz solu- 
tion to the Vlasov equation has been proposed which 
allowed us to explicitly calculate the core-halo distribu- 
tion function for confined plasmas and self-gravitating 
systems |2l, |3| . In this Letter we will show that an ansatz 
solution is also possible for the HMF model. The the- 
ory proposed allows us also to locate the non-equilibrium 
para-to-ferromagnetic phase transition, which earlier the- 
ories incorrectly predicted to be of second order f8]. All 
the results are compared with the molecular dynamics 
simulations performed using symplectic integrator, and 
are found to be in excellent agreement. 



The HMF model consists of N, XY interacting spins, 
whose dynamics is governed by the Hamiltonian 
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where angle 9i is the orientation of the ith spin and 
Pi is its conjugate momentum 's'-llOl- The macroscopic 
behavior of the system is characterized by the magne- 
tization vector M = (MxjMy), where = {cos 9), 
My = (sin 6*), and (• ■ •) stands for the average over all 
particles. The modulus, M = |M| serves as the order pa- 
rameter which measures the coherence of the spin angular 
distribution: for M = we have a completely incoherent 
state, whereas for finite M there is some degree of coher- 
ence. Hamilton's dynamic equations for each spin can be 
expressed in terms of the total magnetization and take a 
particularly simple form 9i — F(6i), where the force on 
each spin is F(^^i) ~ —AIxSm9i + My cos9i. The average 
energy per particle is u — H/N = + (1 — M^)/2. 
Since the Hamiltonian does not explicitly depend on 
time, u is conserved along the temporal evolution. For 
simplicity we will consider initial distributions of "water- 
bag" form in the {9,p) reduced phase-space (^-space). 
Without loss of generality, we choose a frame of refer- 
ence where {9) = and {p) =0. The one-particle initial 
distribution function then reads 
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e(f?o-|0|)e(po- W), 
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where Q is the Heaviside step function, and a-nd 
Pq are the maximum absolute values of the angle and 
the momentum, respectively. Such initial distributions 
are characterized by Mj. ~ Mq, My = 0, and u — 
pI/6 + {1 - M§)/2, where Mq = sin(6'o)/6'o is the ini- 
tial magnetization. Because of the symmetry of /o with 
respect to 9 = 0, in the thermodynamic limit My ~ 
throughout the evolution, so that the macroscopic dy- 
namics is completely determined by Mx{t). 
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As the system evolves, the particle distribution in the 
phase-space changes and eventually reaches a stationary 
state or a limit cycle. If N is finite, the stationary state 
will be described by the equilibrium Maxwell-Boltzmann 
(MB) distribution. In the thermodynamic limit N — oo, 
however, the system becomes trapped in a non-ergodic 
non-mixing state, the life time of which diverges with 
the number of particles. In this limit, the dynamical 
evolution of one-particle distribution function f{9,p, t) is 
governed exactly by the collisionless Boltzmann (Vlasov) 
equation fll'|. 



(3) 



The left-hand side of this equation is just the convec- 
tive derivative of the one particle distribution function. 
Therefore, a collisionless Hamiltonian system evolves 
over the phase space as an incompressible fluid. Fur- 
thermore, the incompressibility implies that during the 
temporal evolution the phase-space density can never ex- 
ceeds the maximum of the initial distribution function. 

Although the MB distribution is also a stationary solu- 
tion of the Vlasov equation, unlike for Boltzmann equa- 
tion, it is not a global attractor of the dynamics, so that 
an arbitrary initial distribution will not evolve to the MB 
equilibrium. The collisionless relaxation described by the 
Vlasov equation for systems with long-range interactions 
is, therefore, much more complex than the collisional re- 
laxation governed by the usual Boltzmann equation for 
systems with short-range forces. 

Vlasov equation is time reversible. Thus, on a fine- 
grained scale, temporal evolution never ends. There is no 
fine-grained attractor for the Vlasov dynamics. In prac- 
tice, however, one can never have infinite resolution, and 
there is a finite maximum precision that one can reach in 
any experiment or a numerical simulation. It is on this 
coarse-grained scale that it appears that the evolution 
has reached a steady state. Unlike for Boltzmann equa- 
tion, however, the stationary coarse-grained distribution 
function depends explicitly on the initial condition. 

Recently it has been observed that for systems with 
long-range interactions, such as self-gravitating clusters 
and plasmas [3, Q , the final stationary state has a pe- 
culiar core-halo structure. The mechanism of core-halo 
formation is very similar to the process of evaporative 
cooling. As the dynamics evolves, macroscopic propa- 
gating density waves are formed. Some particles enter in 
resonance with the macroscopic oscillations gaining large 
amount of energy at the expense of the collective motion. 
This is similar to the mechanism of Landau-damping well 
known in plasma physics IJ]. Resonant particles can 



gain sufficient energy to reach high energy states, thus 
forming a diffuse halo. On the other hand, the loss of 
energy dampens the macroscopic oscillations, so that the 
leftover particles become condensed into the low energy 
states, resulting in a dense core. However, because of the 



incompressibility of the Vlasov dynamics, the core can- 
not completely freeze - i.e., collapse to the minimum of 
the potential energy. Instead, the distribution function 
of the core particles progressively approaches the maxi- 
mum phase space density allowed by the initial distribu- 
tion — all the low energy states become fully occupied 
by the core particles. Although the HMF model appears 
to be very different from either self-gravitating clusters 
or confined plasmas, we find that its dynamical evolution 
follows exactly the same scenario as described above. 

In the case of the HMF, the oscillations of the magne- 
tization M play the role of collective oscillations which 
drive some spins to higher energy states, leading to a 
halo formation. The macroscopic oscillations of M are 
significantly damped in one or two periods of oscillation. 
The extent of the halo is determined on the same time 
scale. As a consequence of the conservation of the total 
energy, the remaining spins populate lower and lower en- 
ergy states, until all of them become fully occupied up 
to the maximum phase space density rjo = l/40oPo- In 
the final stationary state, the core distribution function 
is the same as that of a fully degenerate Fermi gas of 
spin-degeneracy tjq. The core distribution extends up to 
the Fermi energy ep- The value of ej? is yet unknown, 
and must be determined self-consistently. We propose 
an ansatz for the core-halo distribution that describes 
the final (coarse-grained) stationary state reached by the 
HMF model at the end of its dynamical evolution: 

fsie.p) - r;o Mep - e) + xe(e„ - e)e(e - ep)] , (4) 

where e{6,p,Ms) = + 1 — Ms cos 9 is the single- 

spin energy, x is the ratio between the halo and the core 
phase-space densities, Ms is the stationary value of mag- 
netization, and Eh is the maximum energy of the halo 
spins. The energy Eh is determined from the short-time 
dynamics of spins driven by the oscillations of the mag- 
netization. To estimate this value we need an indepen- 
dent equation that describes the dynamical evolution of 
magnetization. We proceed as follows: taking the second 
derivative of M^ we obtain M^ = M^ (sin^ 9) — {p^ cos 9) . 
This equation requires the knowledge of the temporal 
evolution of (sin^ 9) and (p^ cos 9) , which leads to an 
infinite hierarchy of equations. To truncate the hierar- 
chy, we assume that for short times (sin^ 6*) = 1/2 and 
(p^cosS) = (p2)(cos6i) = {2u - 1 + Ml)Mx, where use 
has been made of the conservation of energy, together 
with the condition M = M^. We then find a dynamical 
equation satisfied by the magnetization. 



-M, ( 2u + Ml - - 



(5) 



This equation can be integrated numerically to provide 
the temporal evolution of Mx{t). Since Eq. ([5]) was de- 
rived neglecting the correlations between angles and mo- 
mentums, its validity extends only to very short times. 
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However, the maximum energy of the halo is also deter- 
mined by the very short-time dynamics. Thus, Eq. ([5|), 
should be sufficient to give a reasonable estimate of the 
value of the maximum halo energy. We adopt the follow- 
ing procedure to determine Sh- For a given initial distri- 
bution, we determine the maximum energy attained by 
a group of non-interacting test-spins that are launched 
with the initial conditions selected from the distribution 
function, Eq. ^ . Their dynamical evolution is governed 
by e^ = -M^(t) sm9^ with M^{t) determined by Eq. 
with Ma;(0) = Mq and Mx{t) = 0. We solve this equation 
over a short time corresponding to two periods of oscilla- 
tion of Mx- The Eh-, then, corresponds to the maximum 
energy obtained by any of the test-spins. 

Once Eh has been determined using the test particle 
dynamics, the remaining parameters of the final station- 
ary distribution — ep, Xi ^-^d Mg — are obtained by 
imposing the conditions of conservation of norm and of 
the total energy, as well as the closure equation for mag- 
netization: 

J fs{e,p)d9dp^i, (6) 
j fs{e,p)e{e,p,Ms)dedp = u, (7) 
J f,{e,p) cose dedp = M,. (8) 

The equations above can be analytically evaluated in 
terms of elliptic integrals, forming a closed set of alge- 
braic equations that must be solved numerically to de- 
termine sp, Xi Ei^nd Ms- 

In fig.(Il]) (a and b) we show a snapshot of the phase- 
space and of the distribution of spin energies. The core- 
halo separation can be seen very clearly. The energy Eh 
delimits the particle distribution over the phase space, 
while the Fermi energy restrict the extent of the core 
region. In this example the initial distribution has Mq = 
0.80 and u — 0.55, while the magnetization in the final 
stationary state is Ms = 0.56. In panels (c) and (d) we 
compare the theoretically calculated spin and momentum 
distributions with the molecular dynamics simulations. 
As can be seen, excellent agreement is found between 
the theory and the simulations. 

For all the initial distributions with finite magnetiza- 
tion, we find that if the energy per particle is less than 
u < ui there are two real roots of eq. ([8]), Ms = and 
AIs 7^ 0. The root Ms = is unstable, so that only the 
solution with finite magnetization has a physical mean- 
ing. On the other hand for u > Uu, there is only one real 
root, Ms = 0. For the interval of energies ui < u < Uu 
there are three distinct roots, see fig ([2]). The theory, 
therefore, predicts that there is a first order phase tran- 
sition in the interval between ui < u < Uu- This is 
precisely what is found in simulations, see fig.Q. Un- 
fortunately, differently from the equilibrium phase tran- 
sitions, here we do not have a free energy, which would 




FIG. 1: Snapshots of (a) the phase-space and (b) of spin en- 
ergies £, as a function of angle. The snapshot is taken at 
t = 10000, using N=20000 spins. The solid curves correspond 
to the calculated Fermi energy ef (red) and to the maximum 
halo energy Eh (blue), i.e., e{6,p, Ms) = + 1 — Ms cos 6 = 
Ef and e{9,p, Ms) = Eh, respectively. We see that the Fermi 
energy curve perfectly encloses the high density region. The 
maximum halo energy obtained using the test particle dynam- 
ics and Eq. ((S)) also delimits well the extent of the particle 
distribution in the phase space (blue line). Panels (c) and 
(d) show the angle and the momentum distributions, respec- 
tively. Solid curves are the theoretical predictions obtained 
using the distribution function of Eq. Q, and points are the 
results of molecular dynamics simulations averaged over 20 
runs. The initial distribution has Mq = 0.80 and u — 0.55. 
Note that the present theory predicts that the maximum en- 
ergy attained by any spin will be bounded by Eh, while the- 
ories based on generalized entropies predict that this energy 
distribution is unbounded, d ecay ing either exponentially or 
algebraically 9], but see also [la ]. 



allow us to precisely locate the transition point using 
the Maxwell construction. All we can do is delimit the 
location of the first-order phase transition within a nar- 
row interval ui < u < Uu, shaded gray in fig On 
the paramagnetic side of the phase transition, the sys- 
tems becomes trapped in an out-of-equilibrium limit cy- 
cle associated with appearance of resonance islands in 
the phase space [13], characterized by strong oscillation 
of magnetization around M — 0, see fig. 

We have studied the paradigmatic system with long- 
range interactions the, so called, Hamiltonian Mean-Field 
Model. It is shown that in the thermodynamic limit this 
model does not relax to the usual Maxwell-Boltzmann 
distribution. Instead the final stationary state of the 
HMF has a core-halo distribution function which is an 
ansatz solution to the Vlasov equation. The theory al- 
lows us to explicitly calculate the spin and the momen- 
tum distribution functions, both of which are found to 
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FIG. 2: Magnetization of the final stationary state Ms as a 
function of tlie initial energy per particle u for Mo = 0.40. 
The curves are the theoretical predictions obtained using the 
distribution function of Eq. and the points are the re- 
sults of the molecular dynamics simulations with N = 10® 
spins. The solid curve corresponds to the stable solutions, 
whereas the dashed curve to unstable solutions. The grey area 
corresponds to the metastable region in the parameter space 
where the phase transition must occur. The inset shows the 
result of a large number of molecular dynamics simulations, 
performed in the close vicinity of the phase transition. The 
abrupt change in magnetization, as u is varied indicates a 
first-order phase transition, as is predicted by the theory. For 
all simulations, the integration was performed up to t = 1200. 
The final Ms was obtained by averaging the magnetization 
over an additional time interval At = 800. We note that the 
theory based on Lynden-Bell's coarse-grained entropy incor- 
rectly predicts that the phase transition at this point will be 
of second order 



be in excellent agreement with the simulations, without 
any adjustable parameters. We also find that the HMF 
model undergoes a first-order fcrro-to-para phase transi- 
tion in the region of parameter space where earlier the- 
ories based on the Lynden-Bell (LB) coarse-grained en- 
tropy predicted a second order phase transition [Sjl. It is 
interesting to note that the previous simulations [Sj, |l4 1 
— performed over a fairly short time span t < 150, be- 
fore the system has fully relaxed — failed to notice this 
incorrect prediction of the LB theory. 

The present theory — as well as the earlier results on 
non- neutral plasmas 0, and self-gravitating systems in 



15 1 and two [3| dimensions — suggests that there 



is a significant degree of universality in coUisionless re- 
laxation dynamics. The core-halo distribution function 
appears to be a universal attractor — in a coarse-grained 
sense — of systems with long-range interactions, analo- 
gous to the MB distribution for collisional systems with 
short-range forces. 
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FIG. 3: Dynamical evolution of magnetization on two sides 
of the first-order phase-transition. On the paramagnetic side 
(u = 0.6170), the magnetization oscillates around zero; while 
on the ferromagnetic side (u = 0.6165), the oscillations are 
damped and magnetization converges to Ms predicted by the 
theory. 
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